Library

library(Rfastp)
library(qckitfastq)
library(plotly)
library(plyr)

Data

folder <- "../01_data"
files <- list.files(path = folder, pattern = ".gz", full.names = T)
names(files) <- gsub(pattern = "\\.fastq.gz", replacement = "", x = basename(files))

`Rfastp`

Analysis

qc_SE <- function(file, outdir, name, ad1) {
  res <- rfastp(
    read1 = file,
    outputFastq = paste(outdir, name, sep = "/"),
    thread = 4,
    cutMeanQual = 20,
    cutFrontMeanQual = 20,
    cutTailMeanQual = 20,
    cutLowQualTail = T,
    cutLowQualFront = T,
    qualityFiltering = T,
    maxNfilter = 5,
    qualityFilterPhred = 20,
    qualityFilterPercent = 40,
    lengthFiltering = T,
    minReadLength = 15,
    maxIndexMismatch = 1,
    adapterTrimming = TRUE,
    adapterSequenceRead1 = ad1
  )
  # file.remove(paste0(outdir, "/", name, ".html"))
  return(res)
}

res <- parallel::mclapply(setNames(names(files), names(files)), function(x) {
  qc_SE(file = files[[x]], outdir = "output/", name = x, ad1 = "CTGTAGGCACCATCAAT")
}, mc.preschedule = F, mc.cores = 16)

json_files <- list.files(path = "output", pattern = "json", full.names = T)
names(json_files) <- gsub(pattern = "\\.json", replacement = "", x = basename(json_files))

trim_files <- list.files(path = "output", pattern = "gz", full.names = T)
names(trim_files) <- gsub(pattern = "\\_R1.fastq.gz", replacement = "", x = basename(trim_files))

Reading and getting information from JSON file

General

#' Make general information data frame from
#'
#' @param json A JSON file with QC information
#'
#' @return A data frame with FASTQ file information
#' @export
#'
#' @examples
make_general_df <- function(json) {
  js <- rjson::fromJSON(file = json)

  # pkg <- as.character(packageVersion("Rfastp"))
  seq <- js$read1_before_filtering$total_cycles
  seq_type <- ifelse(test = any(grepl(pattern = "read2", x = names(js))), yes = "paired-end", no = "single-end")
  cycle <- paste0(seq_type, " (", seq, " cycles)")

  mean_len_bf <- js$summary$before_filtering$read1_mean_length
  mean_len_af <- js$summary$after_filtering$read1_mean_length
  dup_rate <- ifelse(test = seq_type == "paired-end",
    yes = paste0(round(x = (js$duplication$rate * 100), digits = 2), "%"),
    no = paste0(
      round(x = (js$duplication$rate * 100), digits = 2),
      "% (may be overestimated since this is SE data)"
    )
  )

  df <- data.frame(
    test = c(
      # "Rfastp version",
      "Sequencing",
      "Mean length before filtering",
      "Mean length after filtering",
      "Duplication rate"
    ),
    value = c(
      # pkg,
      cycle,
      mean_len_bf,
      mean_len_af,
      dup_rate
    )
  )
  return(df)
}

parallel::mclapply(setNames(names(json_files), names(json_files)), function(x) {
  knitr::kable(make_general_df(json = json_files[[x]]), col.names = NULL)
}, mc.preschedule = F, mc.cores = 4)

$SRR955628

Sequencing	single-end (50 cycles)
Mean length before filtering	50
Mean length after filtering	27
Duplication rate	35.69% (may be overestimated since this is SE data)

$SRR955629

Sequencing	single-end (50 cycles)
Mean length before filtering	50
Mean length after filtering	27
Duplication rate	35.39% (may be overestimated since this is SE data)

$SRR955633

Sequencing	single-end (50 cycles)
Mean length before filtering	50
Mean length after filtering	29
Duplication rate	45.74% (may be overestimated since this is SE data)

$SRR955634

Sequencing	single-end (50 cycles)
Mean length before filtering	50
Mean length after filtering	26
Duplication rate	26.43% (may be overestimated since this is SE data)

Before QC

make_before_filt_df <- function(json) {
  js <- rjson::fromJSON(file = json)
  tot_reads <- js$read1_before_filtering$total_reads
  tot_bases <- js$read1_before_filtering$total_bases
  q_20_bases <- js$read1_before_filtering$q20_bases
  q_20_rate <- js$summary$before_filtering$q20_rate * 100
  q_30_bases <- js$read1_before_filtering$q30_bases
  q_30_rate <- js$summary$before_filtering$q30_rate * 100
  gc <- paste0(round(x = js$summary$before_filtering$gc_content * 100, digits = 2), "%")

  df <- data.frame(
    test = c(
      "Total reads",
      "Total bases",
      "Q20 bases",
      "Q30 bases",
      "GC content"
    ),
    value = c(
      tot_reads,
      tot_bases,
      paste0(q_20_bases, " (", q_20_rate, "%)"),
      paste0(q_30_bases, " (", q_30_rate, "%)"),
      gc
    )
  )
  return(df)
}


parallel::mclapply(setNames(names(json_files), names(json_files)), function(x) {
  knitr::kable(make_before_filt_df(json = json_files[[x]]), col.names = NULL)
}, mc.preschedule = F, mc.cores = 4)

$SRR955628

Total reads	21354015
Total bases	1067700750
Q20 bases	1035533864 (96.9873%)
Q30 bases	972278842 (91.0629%)
GC content	49.44%

$SRR955629

Total reads	15145094
Total bases	757254700
Q20 bases	732372019 (96.7141%)
Q30 bases	685909354 (90.5784%)
GC content	49.9%

$SRR955633

Total reads	18771501
Total bases	938575050
Q20 bases	910665276 (97.0264%)
Q30 bases	853893917 (90.9777%)
GC content	49.61%

$SRR955634

Total reads	8999053
Total bases	449952650
Q20 bases	433546624 (96.3538%)
Q30 bases	405844827 (90.1972%)
GC content	52.72%

After QC

make_after_filt_df <- function(json) {
  js <- rjson::fromJSON(file = json)
  tot_reads <- js$read1_after_filtering$total_reads
  tot_bases <- js$read1_after_filtering$total_bases
  q_20_bases <- js$read1_after_filtering$q20_bases
  q_20_rate <- js$summary$after_filtering$q20_rate * 100
  q_30_bases <- js$read1_after_filtering$q30_bases
  q_30_rate <- js$summary$after_filtering$q30_rate * 100
  gc <- paste0(round(x = js$summary$after_filtering$gc_content * 100, digits = 2), "%")

  df <- data.frame(
    test = c(
      "Total reads",
      "Total bases",
      "Q20 bases",
      "Q30 bases",
      "GC content"
    ),
    value = c(
      tot_reads,
      tot_bases,
      paste0(q_20_bases, " (", q_20_rate, "%)"),
      paste0(q_30_bases, " (", q_30_rate, "%)"),
      gc
    )
  )
  return(df)
}


parallel::mclapply(setNames(names(json_files), names(json_files)), function(x) {
  knitr::kable(make_after_filt_df(json = json_files[[x]]), col.names = NULL)
}, mc.preschedule = F, mc.cores = 4)

$SRR955628

Total reads	15818144
Total bases	438765815
Q20 bases	434345957 (98.9927%)
Q30 bases	412966207 (94.12%)
GC content	47.39%

$SRR955629

Total reads	9752048
Total bases	263972256
Q20 bases	261287984 (98.9831%)
Q30 bases	248299061 (94.0626%)
GC content	47.46%

$SRR955633

Total reads	12578898
Total bases	366313296
Q20 bases	361950551 (98.809%)
Q30 bases	342077299 (93.3838%)
GC content	47.02%

$SRR955634

Total reads	1746985
Total bases	45581788
Q20 bases	45018894 (98.7651%)
Q30 bases	42379168 (92.9739%)
GC content	49.5%

Filtering result

make_filt_df <- function(json) {
  js <- rjson::fromJSON(file = json)

  reads <- js$summary$before_filtering$total_reads
  read_pass <- js$filtering_result$passed_filter_reads
  read_pass_pct <- round((read_pass / reads) * 100, 2)
  read_low_Q <- js$filtering_result$low_quality_reads
  read_low_Q_pct <- round((read_low_Q / reads) * 100, 2)
  read_with_N <- js$filtering_result$too_many_N_reads
  read_with_N_pct <- round((read_with_N / reads) * 100, 2)
  read_short <- js$filtering_result$too_short_reads
  read_short_pct <- round((read_short / reads) * 100, 2)
  read_long <- js$filtering_result$too_long_reads
  read_long_pct <- round((read_pass / reads) * 100, 2)

  df <- data.frame(
    test = c(
      "Reads passed filters",
      "Reads with low quality",
      "Reads with too many N",
      "Reads too short",
      "Reads too long"
    ),
    value = c(
      paste0(read_pass, " (", read_pass_pct, "%)"),
      paste0(read_low_Q, " (", read_low_Q_pct, "%)"),
      paste0(read_with_N, " (", read_with_N_pct, "%)"),
      paste0(read_short, " (", read_short_pct, "%)"),
      paste0(read_long, " (", read_long_pct, "%)")
    )
  )
  return(df)
}


parallel::mclapply(setNames(names(json_files), names(json_files)), function(x) {
  knitr::kable(make_filt_df(json = json_files[[x]]), col.names = NULL)
}, mc.preschedule = F, mc.cores = 4)

$SRR955628

Reads passed filters	15818144 (74.08%)
Reads with low quality	12166 (0.06%)
Reads with too many N	0 (0%)
Reads too short	5523705 (25.87%)
Reads too long	0 (74.08%)

$SRR955629

Reads passed filters	9752048 (64.39%)
Reads with low quality	9227 (0.06%)
Reads with too many N	0 (0%)
Reads too short	5383819 (35.55%)
Reads too long	0 (64.39%)

$SRR955633

Reads passed filters	12578898 (67.01%)
Reads with low quality	11657 (0.06%)
Reads with too many N	0 (0%)
Reads too short	6180946 (32.93%)
Reads too long	0 (67.01%)

$SRR955634

Reads passed filters	1746985 (19.41%)
Reads with low quality	3238 (0.04%)
Reads with too many N	0 (0%)
Reads too short	7248830 (80.55%)
Reads too long	0 (19.41%)